function k= zero_eigs_adj(k, p, h)

opts = optimset('tolfun', 1e-10);
k = fsolve(@(x) fun(x, h, p), k, opts);

fprintf('k = %.4e\n', k);

function F = fun(k, h, p)


A = w_coeffs(1, h, k, p);
w = @(z)  (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);

w = @(z) p.Ma * 1e3 * z.^2 .* (h - z);

v_0 = (1 - p.beta) * (1 - 1 / h);
c = p.beta + v_0 - p.delta * (1-p.beta)*(p.beta + v_0) / 2;
gamma = p.delta * (1 - 2 * c);


% the problem in general

F = -k * sinh(k*h) - gamma * cosh(k*h) + p.delta * (1 - p.beta) * (p.beta + v_0) / h^2 * quad(@(z) w(z) .* z .* cosh(k*z), 0, h, 1e-10);